NUMERICAL INVESTIGATION OF THE SMALLEST EIGENVALUES OF 
THE p-LAPLACE OPERATOR ON PLANAR DOMAINS 



JIRl HORAK 

Abstract. The eigenvalue problem for the p-Laplace operator with p > 1 on planar domains 
with the zero Dirichlet boundary condition is considered. The Constrained Descent Method and 
the Constrained Mountain Pass Algorithm are used in the Sobolev space setting to numerically 
investigate the dependence of the two smallest eigenvalues on p. Computations are conducted 
for values of p between 1.1 and 10. Symmetry properties of the second eigenfunction are also 
examined numerically. While for the disk an odd symmetry about the nodal line dividing 
the disk in halves is maintained for all the considered values of p, for rectangles and triangles 
symmetry changes as p varies. Based on the numerical evidence the change of symmetry in this 
case occurs at a certain value po which depends on the domain. 



1. Introduction 

For a bounded domain Q C M. N , N G N and a parameter p G (l,oo) consider the nonlinear 
eigenvalue problem 

—A p u = \\u\ p ~ 2 u in Q, 
u = on dtt 

to be solved for a real function u : Q —> R and a parameter A G R. The operator A p u := 
div (|Vu| p-2 Vii) is called the p-Laplace operator. If for a certain A a nontrivial weak solution 
u G Wq ,p (Q) of (pQ) exists, we call A and u Dirichlet eigenvalue and eigenfunction of the p-Laplace 
operator, respectively. Problem (pQ) is homogeneous but in general not additive in u. 

From [TJ |20j [2TJ |4] and others it is well known that there exists the smallest eigenvalue Ai and 
that it is positive, isolated and simple (i.e., the corresponding eigenfunction u\ is unique up to 
multiplication by a constant). Moreover, for any eigenfunction u it holds: u corresponds to Ai 
if and only if it does not change its sign on Q. In jT2] using a variational approach the authors 
constructed a nondecreasing sequence of eigenvalues accumulating at infinity. Since between Ai 
and the next member of this sequence there are no other eigenvalues, as it was shown in [2j, we call 
this second smallest eigenvalue A2 and a corresponding eigenfunction i£ 2 . In general, however, it 
is not known yet whether this sequence contains all the eigenvalues. Nodal domains of variational 
eigenfunctions were studied in The regularity results of [10J imply that any eigenfunction 

(perhaps redefined on a class of measure zero) is of class C 1,a (Q) for some a > 0. 

An early attempt at computing several eigenpairs of the p-Laplace operator on a planar do- 
main (TV = 2) numerically is due to Brown and Reichel |7J. Under the assumption of radial 
symmetry they used a shooting method for the resulting ordinary differential equation. The first 
genuinely two-dimensional approach was taken by Yao and Zhou [24J using their local minimax 
method based on a variational formulation. For a square £1 = {(^i,^) | #i,#2 G (0,2)} and 
p G {1.75,2.5,3.0} the authors computed approximations to seven eigenvalues and corresponding 
eigenfunctions. They observed that the found eigenfunction U2 has an odd symmetry about x\ = 1 
for p < 2 and about x\ = £2 for p > 2. 

The goal of the current work is to apply the numerical variational methods of [14J to compute 
approximations of the two smallest eigenvalues and to visualize the corresponding eigenfunctions 
on a planar domain. In particular the focus is 
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• to extend the Constrained Mountain Pass Algorithm from the Hilbert space setting (as 
described in [8] and [14]) to the Banach space W jP (Q); to verify that this algorithm is 
suitable even for computations with p "far" from 2; 

• to observe the behavior of the eigenpairs for a large range of p and compare it with the 
known theoretical results about the asymptotics for p — >> 1 and p — > oo; 

• to observe changes in symmetry of u<i on various domains. 

In Section [2] we review known results about the variational properties of Ai and A2 and their 
asymptotic behavior. The variational numerical methods applied to compute the eigenpairs are 
summarized in Sec. [5] The choice of a descent direction in the Banach space W^ p (yt) is discussed 
in detail here, too. In Sec. [4] we present the numerical results for several planar domains. We 
pay a particular attention to the dependence of the eigenvalues on p and changes of symmetry of 
the second eigenfunction. Several issues concerning the application of the numerical methods (like 
mesh refinement, choice of parameters, etc.) are addressed in Sec. [5] Finally, Sec. [6] summarizes 
our numerical observations and the Appendix provides proofs of several claims used in Sec. [3j 

2. Background material 

In the Introduction we mentioned the existence of the first two eigenvalues Ai and A2. Now we 
review some known results about their variational characterization based on the above references 
and their asymptotic behavior for p close to 1 and p large. 

2.1. Variational characterization of Ai and A2. Define two continuously Frechet differentiable 
functionals J, J G C 1 (W^ p (Q) ,R) : 

(2) I(u) := / \Vu\p dx, J(u) := f \u\ p dx. 

Jq Jq 

Their Frechet derivatives I'(u), J'{u) are members of the dual space of W c \' p (tt) which we denote 
by W~ 1,q (Q), where ^ + ^ = 1, and are given by 

(3) = p [ \Vu\ p - 2 VuV(j) dx, (J'(u)^) = p f \u\ p - 2 u(j) dx. 

Jq Jq 

Two observation can be made: 1. After testing (pQ) with G Wq ,p (Q) and integrating by parts 
it becomes clear that (pQ) is the Euler-Lagrange equation T \u) — \J'(u) = (up to the factor p) 
which all critical points of / with respect to the constraint 

(4) S := {u e w^ p {n) I J(u) = l} 

must satisfy for some value of the Lagrange multiplier A. 

2. If (\,u) is an eigenpair and we test ([!} with u, we obtain the Rayleigh quotient 

(5) A= J>^! 

In M P dx ' 

Since both its numerator and denominator are homogeneous of the same degree in u, finding the 
smallest eigenvalue A is the same as minimizing / on S: 

(6) Ai = minJ(^i). 

A variational minimax characterization of the second eigenvalue A2 based on the Krasnoselskii 
genus was given in [12]. Alternatively, since for u\ G S both u\ and — u\ are local minimizers of / 
on S, a mountain pass characterization of A2 is also possible |9J: 

(7) A2 = inf max I(u). 

7^rnG7([0,l]) 

where T = {7 G C([0, 1],S) \ 7(0) = ui, 7(1) = — ui} is the family of all paths in S connecting 
the two local minimizers. Hence for the numerical computations we have the setting required by 
the Constrained Mountain Pass Algorithm of |14| . 
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2.2. Asymptotic behavior of Ai and A2 as p — >> 1. To make the dependence of an eigenvalue 
on the domain Q and the parameter p explicit in our notation we will write \(Q;p) if necessary 
(and similarly for eigenfunctions) . The main result of [T7] implies that for Q with a Lipschitz 
boundary 

Per(D) 

(8) lim Ai(0;p) = MO), where fti(fi) := min . }. J 
p^l ' ' DdVt \D\ 

is called Cheeger constant, Per(D) denotes the perimeter of D measured with respect to W N and 
\D\ its TV-dimensional Lebesgue measure. A minimizer in the definition of h\(Q) is called a Cheeger 
set of ft. Furthermore, any convex planar domain possesses a unique Cheeger set Cq and 

(9) lim ui(Q;p) = xc& m L l along a subsequence. 

Here the eigenfunctions ui(Q]p) have been normalized to 1 in the L°°-norm, xc n is the indicator 
function of Cq. 

A detailed description of how to find the Cheeger set Cq for a convex planar domain Q is given 
in |18| . Its main property is 

(10) Cq = [j j 5 C tt B is a ball of radius | . 
In [22J it was shown that for Q with a Lipschitz boundary it holds: 

(11) limA 2 (ft;p) = /i 2 (ft), 

p—ti 



where 

(12) h 2 (tt) := min 1/iG 



3D U D 2 C O, D x H D 2 = and max < /i 

1=1,2 IDA 



is called the second Cheeger constant and the convention Per(D)/\D\ — 00 is used if \D\ = 0. Any 
two sets Di,D 2 for which the minimum in the definition of h 2 {Cl) is achieved are called coupled 
Cheeger sets of ft. For a result about the L 1 -convergence of the second eigenfunctions we refer to 
[22| Thm. 5.11]. 

2.3. Asymptotic behavior of Ai and X 2 as p — >• 00. For a bounded domain ^ of M N a limit 
problem of (pQ) as p — >• 00 is studied in [T6| [15] for an unknown function u and an unknown real 
parameter A (see [15l Definition 2.1]). The smallest A for which this limit problem admits a 
nontrivial viscosity solution is called the first oo-eigenvalue and denoted Ai. For Ai there exists a 
positive viscosity solution and it holds: 

(13) lim (X 1 (n;p)) 1/P =Ai(n), 

(14) Ai(^) = — , where T\ := sup{r > | 3 an open ball B C Vt of radius r}, 

(15) A 1 (0)=min{ 11 ^ 11 ^ „ € \ {0}} . 

The characterization ([15]) is an analogy of ([5]) and ©. Furthermore, for any sequence {^(fJ;^)}?^ 
with pi —> oo and ||^i(^;Pi)||LPi(n) = 1 there exists a subsequence converging uniformly to a vis- 
cosity solution of the limit problem for Ai(Q). 

The smallest A for which the limit problem admits a viscosity solution with at least two nodal 
domains is called the second oo-eigenvalue and denoted A2. From the definition it follows that 
Ai < A2. If Ai < A2, then for A G (Ai, A2) zero is the only solution of the limit problem. It holds: 

(16) lim (X 2 (n;p)) 1/P = A 2 (a), 
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(17) A2(^) = — , where r2 := sup{r > | 3 disjoint open balls B\, B2 C ft of radius r}, 

(18) A 2 (ft) = inf max ||Vtx|| L ao (n) , 

7^r n G7([0,l]) 

where T is defined as in ([7J), ui is any first oo-eigenfunction and S := {u £ W lCC (Q) | |M|L°°(n) = 
1}. Furthermore, for any sequence {^(ft;^)}^ with pi — >> 00 and ||i£2(ft;Pi)IU p *( fi ) = 1 there 
exists a subsequence converging uniformly to a viscosity solution of the limit problem for A 2 (ft) 
which has at least two nodal domains. 

3. Numerical methods 

An overview of the numerical methods used to compute approximations of the first and the 
second Dirichlet eigenpair of the p-Laplace operator is given in Fig. [T] In this section we will 
describe these methods. Our goal is to find u\ as a minimizer of I on S according to (|SJ) and 
U2 as a mountain pass point of I on S according to ([7j). We first discretize the planar domain 
ft using a mesh of triangles and apply the finite element method to approximate Wq ,p (Q) by a 
finite dimensional subspace. Then we fix p G (1, 00) and use a variant of the Constrained Steepest 
Descent Method (CDM) to find the first eigenpair, and the Constrained Mountain Pass Algorithm 
(CMPA) to find the second eigenpair. We implement both methods based on [14J. There are, 
however, several important issues arising from the fact that we work in a Banach space and not a 
Hilbert space as in |14| . How to deal with these issues will also be explained in this section. For 
the computation of the descent direction the Augmented Lagrangian Method of [13] is applied. 



loop over a range of values of p G (1, 00) 
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Figure 1. Flowchart of the numerical computations. 



3.1. Finite element method. A finite element approximation of the p-Laplacian was studied in 
[3j. We adopt this approach for our computations. The planar domain ft is approximated by a 
polygonal domain Q h which is partitioned into a finite number of triangles of diameter at most h. 
Let {di}i =1 be the set of those triangle vertices which lie in the interior of Q h . Functions {</>i}i =1 
forming a basis of the /c-dimensional subspace Vq of W^ p {£l h ) are chosen linear on each triangle 
with (f>i(aj) = 5ij, where Sij is the Kronecker delta, and zero on dQ h . The space Vq is our finite 
element approximation of the Sobolev space Wq' p (Q). 

In |3J a detailed description of this method was given for the boundary value problem 

-A p u = f in ft, 
^ ' u = on dn 

with the right-hand side / G L 2 (Q). Since it is a straightforward task to adapt it to our problem 
(PQ) with A|^| p_2 ^ on the right-hand side we will not show the details here. 

We will however mention one additional technical detail involved. The evaluation of the func- 
tional given in ([2]) and ([3j) for functions from Vq amounts to adding up the contributions of the 
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individual triangles that make up Q h . For example, for I(u) with u G Vq one merely needs to 
integrate a constant on every triangle. The situation is different for J(u). Let T C Q h be a triangle 
with area |T| and vertices A, B, and C and let u be a linear function on T with values ua, ub, 
and at these vertices, respectively. If these values are mutually different, then the following 
formula holds: 

(20) f\uf dx 2 ' r ' (\uc\^-\us\^ Kr 2 -M— 



(p + l)(p + 2)(uc - ma) V u c -u B u A - u B 

By inspecting this formula we see that great care must be taken when implementing it to avoid 
numerical cancellations. This is crucial for the success of our method. A similar situation occurs 
when evaluating (J f (u),(j)) (or the right-hand side of equation in (pQ) in the weak formulation). 

3.2. Direction of descent. An important ingredient of the variational numerical methods CDM 
and CMPA is finding a descent direction of the functional / on the constraint set S. How this is 
accomplished in the Hilbert space setting was shown in [14J: Let \7I(u) be the Riesz representation 
of I'(u) (i.e., the gradient) and P u the orthogonal projection on the tangent space of S at u G S. 
Then 

(21) w u = -P u VI{u), ueS 

gives the steepest descent direction of / at u with respect to S. 

Because of the lack of orthogonality in the Banach space Wq ,p (Q) we need to take a different 
approach. Let 

(22) T U S :={ve W^ P (Q) | (J'(u),v) = o} , |MI := Qf |V^ X 

denote the tangent space of S at u G S and the norm of v G Wq' p (^), respectively. The problem 
of finding the steepest descent direction of / with respect to S can be written as follows: for a 
given u G S which is not a critical point of / with respect to S 

(23) minimize (I'(u),w) subject to w G {v G T U S \ \\v\\ = 1}. 



It has a unique solution as Lemma [A. 1| in the Appendix shows. The Euler-Lagrange equation that 
this solution must satisfy can be written in the form 

(24) - A p w = p (-A p u - a\u\ p - 2 u) , 

where a,/3 G R are unknown. The coefficient j3 comes from the requirement ||u>|| = 1. After 
testing ([24]) by w it can be seen that j3 < since the minimum in ([23]) is negative. For p ^ 2 
finding the right a is not an easy problem. 

We will try to find a different convenient descent direction instead, not necessarily the steepest 
one. A simple calculation shows that under no constraints the steepest descent direction of I at 
u G B is given by —u. For u G S we consider the point w u G T U S closest to —u, i.e., the unique 
solution of the minimization problem 

(25) minimize \\w + u\\ subject to w G T U S. 
The minimizer must satisfy the Euler-Lagrange equation 

(26) - A p (w + u) = a\u\ p ~ 2 u 

for some aGi Unlike ([24]) . this equation can be solved easily for w: 

(27) w u = -u + Uulp ! 2uvudx v u , where := (-A,)" 1 {H^u) . 



The operator (— A p ) _1 is discussed later in Sec. 3.3 Lemma A. 2 in the Appendix shows that 
w u is, indeed, a descent direction of I with respect to S. This descent direction is used in our 
implementation of the variational numerical methods CDM and CMPA. 
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Remark. 1. Observe that if — A p were linear, equations ([24]) and ([26]) would coincide (after setting 
/? = —1). Hence in case p = 2 they yield the same descent direction (which is the one given by 

(ED)). 

2. With /3 = — 1 in ([24]) , both equations ([24]) and ([26]) yield a zero solution if and only if u is a 
critical point of / with respect to S. 

3.3. Inverse of the p-Laplace operator. A classical result (see, e.g., [23] Theorem 1.3]) says 
that for any / G W-^ q (Q), the dual of W^ p (n) with ± + ± = 1, the problem 

—A p u = f in ft, 
( 28 ^ u = on ffl 

has a unique weak solution in Wq' p (Q). This means that the operator 

(29) - A p : W^ p (n) W~ hq (Q) given by (-A p u,v) = [ \ Vu\ p ~ 2 VuVv dx 

is invert ible. We denote its inverse by (— A p ) _1 . 

In order to use the descent direction given by ([27]) we need to compute v u first, i.e., we need 
to solve problem ([28]) numerically. For that we apply the Augmented Lagrangian Method of |13| . 
Here we give a brief description of this method. Let V h be again the subspace of continuous 
functions of W^ p (Q h ) which are linear on every triangle of a triangulation of Q h , D h the space 
of functions with values in R 2 defined on Q h which are constant on each triangle, and r > a 
parameter. For the Augmented Lagrangian 

(30) C r (v,t,fjL) = - I \t\ p dx - {f,v) + V - I \Vv-t\ 2 dx + [ p - (Vv - t) dx, 

P Jn 2 J n J n 

where v G V h and t, /i G D h , a saddle point (u, s, rf) is searched for such that 

(31) C r (u, S, /i) < C r (u, S, T]) < C r (v, t, T]) V(v, t, /i) G Vq 1 X X 

A sequence (i/ n ), s( n ), r/ n )) approximating (u,s,rj) is constructed as follows: choose (s^ ),^ 1 )) G 
D' 1 x D h and for n G N solve 

-r A^ (n) = / + V M • V - r s^" 1 ) -V in ft, 

( 32 ) n 

w (n) =0 on 5ft, 

(33) |s (n) | P "V n) +rs (n) = rW n) +r/ n) , 

(n+l) = ^(n) _,_ r (\j lt {n) _ g {n)\ 



(34) r/ n+1 ) =ryW+r(Vtx 



For given s( n_1 ) and the boundary value problem ([32]) can be solved for u^ n \ For this one 
just needs some standard algorithm for finding the inverse of the Laplace operator with Dirichlet 
boundary conditions. The equation in ([32]) is understood in the weak sense: for example the term 
evaluated as J Q rf n ^ ■ V<j)dx for a test function <p G Vq 1 . 

Next, equation ([33]) is used to find s^ n \ The R 2 -norm of must satisfy 

(35) |s (n) | P-1 +r|s (n) | = | r y^ n ) -f r/ n) |, 

which on each triangle is just a scalar nonlinear equation with one unknown. For each triangle it 
can be solved, e.g., by Newton's method. After |s( n )| has been obtained, can be computed 
immediately from ([33]) . 

At the end 77 is updated according to ([34]) and a new iteration step can be started. 

The convergence of this method was studied in [T3] . We use the norm of Vu^ —s^ to measure 
the convergence. 
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(36) —u(t)=w u (t), u(0) = e eS, 



3.4. Constrained Descent Method. The Constrained Descent Method (CDM) is applied to 
find the first eigenpair of the p-Laplace operator: u\ is found as the minimizer of / with respect 
to S, Ai = I(ui). As mentioned above, it differs from the Constrained Steepest Descent Methods 
of [14] in the way the descent direction is chosen. 

The method solves numerically the following initial value problem: 

d 

where w u is given by ([27]) for u G S. Proposition |A.3| in the Appendix states that this problem 
has a unique solution u(t) G S for t G (0, oo) and that u(t) gets arbitrarily close to a critical point 
of / with respect to S as t —> oo. 

After choosing the starting point eo G S and setting := eo the initial value problem 
is solved by repeating the following two steps: First (Euler's step), given i^ n_1 ) find = 
^( n_1 ) -f At^ w u (n-i) with some small value At^ > 0. Second (scaling), define = cu^ n \ 
where the coefficient c G R is chosen such that G S. In case /(^ n )) > I(u( n -^), halve 
the step At^ and compute and again. If this halving has to be repeated and At^ 
becomes very small (smaller than a prescribed threshold value), stop the algorithm. The norm 
of the descent direction ||iu n (™-i) || is used to measure convergence of to an eigenfunction u. 
When computing w u according to ([27)) the integral v := J Q \u\ p ~ 2 uv u dx has to be evaluated. If 

||w u (n-i) || is small, then (l/z/ n_1 )) P 1 approximates the eigenvalue. 

We note that at every step of CDM the Augmented Lagrangian Method of Sec. |3.3| has to be 
applied to compute the descent direction w u ( n -i). 

3.5. Constrained Mountain Pass Algorithm. Suppose that an approximation of the first 
eigenvalue Ai and eigenfunction u\ of the p-Laplace operator have been computed. Constrained 
Mountain Pass Algorithm (CMPA) is applied to find the second eigenpair: is found as a 
mountain pass point of I on S lying "between" the two local minimizers u\ and — u\, X2 = I(v>2). 
Again, it differs from CMPA described in detail in [14J in the choice of the descent direction. 

We give a short summary here based on the original description of the Mountain Pass Algorithm 
by Choi and McKenna [8]: Take a discretized path {zj}jLo c $ connecting zq := ui with zp := 
—u\. After finding the path point z m =: z m& * at which / is maximal along the path, move this 
point a small distance in the tangent space to S at ^ max in the descent direction w z m ax and then 
scale it (as in CDM) to come back to S. Thus the path has been deformed on S and the maximum 
of I lowered. Repeat this deforming of the path until the maximum along the path cannot be 
lowered anymore: a mountain pass point of / with respect to S has been reached. 

To construct the initial path connecting u\ and — u\ in S we choose an intermediate point 
6m G S \ {±^1}, set k := [P/2] and define: 
j 

Zj := u x + -(e M - ^1) for j G {0, . . . , k}, 

Zj := e M + p_ ~ e M) for j G {k,. .. ,P}, 



Zj := CjZj G S (scaling to S as in Sec. 3.4) for j G {0, . . . , P}. 

Connecting u\ and — u\ by a line segment without the intermediate point eM would not work. 
Such a line segment passes through and hence cannot be scaled to get to S. 

Finally, as in CDM, ||i^max|| is used to measure convergence to an eigenfunction u. The cor- 
responding eigenvalue A is computed as in CDM, too. At every step of CMPA the Augmented 



Lagrangian Method of Sec. |3.3| has to be applied to compute the descent direction w z m^. 

4. Numerical results 

In this section numerical results will be given for the following planar domains: the unit disk, 
the square with side length 2, the rectangle with sides 2 and 7/4, the isosceles triangle with base 
and height 1, the isosceles triangle with base 1 and height 3/4, and the equilateral triangle with 
side 1. Unless explicitly stated otherwise the computed eigenfunctions will be plotted as a surface 
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over the domain with heights given by the function values and as a contour plot of these values 
(like, e.g., in Fig. [2|. In order to better compare the shapes the eigenfunctions in these figures 
have been scaled to have the same maximum value. We do not explicitly differentiate between 
two eigenfunctions u and u if u{x) = cu(Tx), where c G M is a scaling coefficient and T : Q — >> ft is 
some symmetry transformation of ft (e.g., for a square a rotation by it/ 2 about the center of the 
square). 

4.1. Unit Disk. Let 

Q = {(x u x 2 ) eR 2 \xl + xl <l} . 

Before presenting the numerical results we make a remark about the radially symmetric case. It 
is known that the first eigenfunction for the disk is radially symmetric. One important question 
about the second eigenfunction for the disk has been whether it is radially symmetric, too. In 
[221 [5] the authors proved that for p close to 1 the answer is no. The eigenvalue problem (pQ) under 
the assumption of radial symmetry u = u(r), r G (0, 1) becomes 

- (r\u f \ p - 2 u f )' = \r\u\ p - 2 u, 
37 V 1 1 J 11 

i/(0) = 0, u(l)=0. 

This and a related problem are treated, for example, in |7J and [5 J, where numerical approaches 
play an important role. For our numerical investigation we adapt the genuine 2D method of Sec. [3] 
in the following ways: 

• all integrals are one-dimensional, 

• the weight r is introduced, 

• the natural boundary condition is implemented at r = (the zero boundary condition 
stays at r = 1). 

Since these modifications are rather elementary, we will not describe them in more detail. We will 
refer to this method as radial 2D method. 

For the computations carried out by the genuine 2D method the domain Q was approximated 
by a polygon and discretized using 68,608 triangles. For the computations carried out by the 
radial 2D method the interval (0, 1) was divided into 1,000 subintervals of the same length. 

Figures [2] and [3] show the eigenfunctions u\ and u 2 computed by the genuine 2D method for 
several values of respectively. The corresponding eigenvalues Ai and A2 for these and other 
values of p are listed in Table [TJa). Figure [5|a) shows the shape of the intermediate point eM on 
the initial path connecting u\ and — u\ we used for CMPA to find u 2 for all the listed values of 
p. The function u 2 found this way seems to posses an odd symmetry with respect to its nodal 
line. The slope of this nodal line in the coordinate system {xi,x 2 ) depends on the computation. 
For the depiction in Fig. [3] we rotated Q in each case to make the slope appear the same. CMPA 
needed between 120 and 600 iterations to converge. 

Figure [5|b) shows an alternative shape of eM- With such an initial path CMPA converged for 
p = 1.1 and p = 1.2 to a radially symmetric function we call u T 2 d (but for higher values of p to 
the oddly symmetric function u 2 ). Figure [IJa) shows u r 2 ad for p = 1.1. 

Figure [5^c) shows yet another choice of eM (radially symmetric). With this intermediate point 
of the initial path and for p = 1.3 and p = 1.4 (but not larger) CMPA seems to converge to a 
radially symmetric function first but after many iterations the path slips down and the algorithm 
converges eventually to the oddly symmetric u 2 . The graph in Fig.[5|d) shows how the maximum 
value of the Dirichlet functional / along the path develops during the run of the algorithm (for 
p = 1.3). The horizontal axis shows the number of iterations. The flat part between iterations 70 
and 260 indicates that the path is staying close to a critical point. When now the norm of the 
descent direction w zmax given in ([27]) computed at the "highest" point £ max of the path gets small 
enough, we stop the algorithm and save this highest point. Since it displays a radial symmetry, 
we call it u 2 ad again. 

The eigenvalues A^ corresponding to the found u 2 ad are also listed in Table [lja). 
Figure Qb) shows profiles of the eigenfunction u 2 ad computed by the radial 2D method for 
several values of p. The eigenvalues Ai and A^ computed by this method for these and other 
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Figure 2. The numerically computed first eigenfunction u\ for the disk. 




p = 1.1 p = 1.4 p = 2.0 p = 2.5 p = 8.0 

Figure 3. The numerically computed second eigenfunction U2 for the disk. 





p = 2.5 



p = 1.1 










p= 8.0 











1/3 1/2 2/3 



Figure 4. The numerically computed radially symmetric second eigenfunction 
u r 2 ad : (a) using the genuine 2D method; (b) using the radial 2D method. The 
profile of u r 2 ad for the radial coordinate r G (0, 1) is shown, the scaling along the 
vertical axis is chosen such that |m ad ||p = 1. 
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p 


Ai 


A 2 


\ rad 


(b) 


P 


Ai 


\ rad 
A 2 


1.1 


2.5690 


4.2008 


5.6809 




1.1 


2.5688 


5.6762 


1.2 


2.9654 


5.0707 


7.2277 




1.2 


2.9653 


7.2251 


1.3 


3.3263 


5.9604 


8.9302 




1.3 


3.3260 


8.9279 


1.4 


3.6740 


6.9072 


10.861 




1.4 


3.6739 


10.858 


1.5 


4.0179 


7.9310 




1.5 


4.0177 


13.073 


1.6 


4.3623 


9.0465 




1.6 


4.3621 


15.626 


1.7 


4.7097 


10.266 




1.7 


4.7095 


18.574 


1.8 


5.0618 


11.604 




1.8 


5.0616 


21.982 


1.9 


5.4194 


13.072 




1.9 


5.4192 


25.921 


2.0 


5.7834 


14.683 




2.0 


5.7831 


30.471 


2.1 


6.1542 


16.452 




2.1 


6.1539 


35.725 


2.2 


6.5320 


18.395 




2.2 


6.5317 


41.788 


2.3 


6.9173 


20.527 




2.3 


6.9169 


48.780 


2.4 


7.3102 


22.866 




2.4 


7.3097 


56.836 


2.5 


7.7107 


25.432 




2.5 


7.7102 


66.112 


3.0 


9.8323 


42.460 




3.0 


9.8314 


137.93 


4.0 


14.683 


110.71 




4.0 


14.681 


559.02 


5.0 


20.351 


273.00 




5.0 


20.347 


2,132.7 


6.0 


26.832 


649.47 




6.0 


26.823 


7,822.6 


8.0 


42.210 


3,430.1 




8.0 


42.182 


97,462 


10.0 


60.784 


17,071 




10.0 


60.715 


1.1359- 10° 



Table 1. Eigenvalues for the disk computed numerically by: (a) the genuine 2D 
method, (b) the radial 2D method. 




Figure 5. (a)-(c) Intermediate point cm of the initial path used in CMPA. (d) 
Maximum value of the Dirichlet functional / along the path during the run of 
CMPA with cm shown in figure (c) for p = 1.3. 



values of p are listed in Table [TJb). The convergence of CMPA does not seem to be sensitive to 
the choice of eM in this case. 

By comparing the values of Ai and A^ in Table [TJa) with those in Table |TJb) which were 
computed by the two different numerical methods we observe that their first three digits coincide 
in almost all the cases. Also, the profiles of ui, u T 2 ad are very close for both methods, respectively 
(cf. Fig |4^a) and the top left graph in (b) for and p = 1.1). We conclude that these are 
numerical approximations of the same eigenvalue-eigenfunction pairs. 
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A A 




Figure 6. Dependence of the numerically computed eigenvalues for the disk on 
p. The three cross symbols in the graph on the right mark the values hi(Q) = 2, 
fc 2 (ft) ~ 3.1543, and h r 2 ad (n) = 4. 




! - Ai = 1 

1 1 , 1 , 1 , 1 , p 

1 23456789 10 



Figure 7. Dependence of the numerically computed eigenvalues for the disk 
raised to 1/p on p. 

The behavior of CMPA suggests that although u r 2 ad is a constrained mountain pass point of I 
among radially symmetric functions, it is not a constrained mountain pass point with no assump- 
tion on the symmetry (cf. Fig. [5jd)) . The case of p = 1.1 and p = 1.2 when CMPA with eM from 
Fig. [5^b) converged to a radially symmetric function and the path did not slip off to asymmetric 
functions with lower values of I seems to contradict this. However, we assume that this was caused 
by the "flat" shape of the landscape of / close to u r 2 ad for p close to 1 and by numerical inaccuracies. 

The dependence of Ai, A2, and A^ on p is presented in Figs. [6] and [7] First of all we observe 
that for all the values of p considered the inequality A2 < A^ holds. Hence this is a numerical 
evidence that the second eigenfunction for the disk is not radially symmetric not only for small p 
but for a large range of p. 

Second, we can observe the following asymptotic behavior: 





Ai 


A 2 


A r 2 ad 


lim p ^i + A 


2 


3.1543 


1 


lim^oo AVp 


1 


2 


3 



Theoretical results for Ai and A 2 were summarize in Sec. [2] The values fti(fi), Ai(fi), and A 2 (fi) 
for the disk are easy to compute. In [22J it was proved that /i 2 (0) for the disk equals the first 
Cheeger constant for the half-disk which is approximately 3.1543. We can observe (Fig. [2]) that 
u\ converges to 1 for p — >> 1 as explained in [17] and to the distance function to the boundary for 
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p — » oo as explained in [16]. In Fig. [3] we observe that for p — >• 1 the function is getting close 
to the indicator function of the Cheeger set for the half-disk on each nodal domain. 

In |5J a numerical evidence is given leading to the conjecture for the asymptotic behavior of 
A^ ad given in the above table. Our numerical results (at least for p —> 1) support this conjecture. 
To motivate these values and the profiles of u 2 ad in Fig. [4] we make the following two remarks: 

Remark. 1. For r G (0, 1) let D(r) be the disk of radius r centered at the origin, A(r) =Q\D(r) 
an annulus. It is easy to show that for r = 1/2 both D and A have the same Cheeger constant 
h r 2 ad (Q) := hi(£)(l/2)) = h 1 (A(l/2)) = 4 (see, e.g., [19J for a result about the Cheeger constant of 
an annulus). The function u v 2 d with its profile shown in Fig. [4] seems to get close to the indicator 
function of D(l/2) and A(l/2) on each nodal domain for p — >• 1. 

2. Under the assumption of radial symmetry two largest disjoint disks of the same radius 
inscribed in Q have radius 1/3. Hence we define = jj^ = 3. The function u T 2 d with its profile 
shown in Fig. [4] seems to get close on each nodal domain to a multiple of the function giving the 
distance to the boundary on D(l/3) and A (1/3) for large p. 

4.2. Square. Let 

SI = {(x u x 2 ) eR 2 \x u x 2 G (0,2)} . 

This domain was discretized using 83,968 triangles. Figures [8] and [9] show the eigenfunctions u\ 
and U2 computed for several values of p, respectively. Table [2] lists the corresponding values of Ai 
and A2. 

Various choices of the intermediate path point eM were used to compute u 2 . Only in case p = 2 
different choices of e^i caused CMPA to converged to different functions 112. Since for the square 
the eigenspace corresponding to the second eigenvalue of the Laplace operator is two-dimensional, 
CMPA converges to some member of this eigenspace depending on the shape of the initial path. 
Figure [9] shows one such eigenfunction. However, even for p = 2 this has no influence on the 
computed value of A2. 

We say that a function has symmetry Si (odd symmetry about xi = 1 and even symmetry 
about X2 = 1) if it belongs to 

Si := {u : Q — >• R | u(xi, X2) = —u(2 — xi, #2), u(xi, X2) = u(xi, 2 — #2)}, 

and symmetry S2 (odd symmetry about xi = X2 and even symmetry about xi — 2 — X2) if it 
belongs to 

£2 := {u : Q — >• R | u(xi, X2) = — u(x2, ^i), ^(^i, X2) = ^(2 — ^2, 2 — #i)}. 

As it was observed in [24J, U2 changes its symmetry at p = 2 from Si for p < 2 to £2 for p > 2. Let 
denote the smallest eigenvalue with an eigenfunction belonging to Si where i G {1,2}. The 



p 


Ai 


A 2 


A<s 2 


1.1 


2.3649 


3.7586 


3.8702 


1.2 


2.6934 


4.5012 


4.6179 


1.3 


2.9986 


5.2500 


5.3715 


1.4 


3.2834 


6.0385 


6.1621 


1.5 


3.5611 


6.8835 


7.0053 


1.6 


3.8356 


7.7971 


7.9118 


1.7 


4.1092 


8.7897 


8.8903 


1.8 


4.3830 


9.8708 


9.9490 


1.9 


4.6581 


11.050 


11.095 


2.0 


4.9349 


12.338 


12.338 



p 


Ai 


A 2 


A<Si 


2.0 


4.9349 


12.338 


12.338 


2.1 


5.2139 


13.684 


13.744 


2.2 


5.4952 


15.144 


15.282 


2.3 


5.7791 


16.725 


16.961 


2.4 


6.0658 


18.438 


18.797 


2.5 


6.3552 


20.293 


20.802 


3.0 


7.8452 


32.107 


33.956 


4.0 


11.038 


74.757 


85.447 


5.0 


14.497 


163.59 


205.08 


6.0 


18.194 


343.77 


477.60 


8.0 


26.221 


1,402.1 


2,443.4 


10.0 


34.990 


5,339.0 


11,888 



Table 2. Eigenvalues for the square. 
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p = 1.1 



p = 1.4 



p = 2.0 



p = 2.5 



p = 8.0 



Figure 8. The numerically computed first eigenfunction u\ for the square. 







p = 1.1 



p = 1.4 



p = 2.0 



p = 2.5 



p = 8.0 



Figure 9. The numerically computed second eigenfunction u 2 for the square. 



values of \s i can be computed using CDM on fl with additional boundary conditions u(l, x 2 ) = 
for x 2 G (0, 2) or x) = for x G (0, 2), respectively, or as the first eigenvalue on the half-domain 

Qhalf where 



^ alf := {(x u x 2 ) G R 2 G (0,1), x 2 G (0,2)} , 
{(x u x 2 ) eR 2 \x ± G (0,2), x 2 G (0,xi)}. 



^^ alf := 



Our numerical observations regarding these eigenvalues and A2 are summarized in Table [3ja) 
and the computed values are listed in Table [2] We stress that A2 was computed with no a priori 
assumptions on symmetry. The dependence of the eigenvalues Ai, A2, Xs ± and A,s 2 on p is further 



plotted in Figures [10] and [TT| 

These figures and Table |3[b) also explain the asymptotic behavior as p — > 1 and p — » 00. While 
the table shows the limit values as given by the theory, the graphs indicate convergence to these 
values (at least for p — >> 1; for p — >> 00 it seems a larger range of p would be needed). The Cheeger 
constants hi shown in the first row of Table [3^b) have been computed according to [17J, [18J by 

(38) fti((0,a) x (0,6)) = — 7====^ for a, 6 > 0. 



a + b - y/(a - hf 



nab 



The evaluation of Ai in the second row is straightforward. 
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A<S 2 


(b) 




Ai 


A^i 


A<s 2 


p < 2 


= A 2 


> A 2 




lim p ^i + A 




4-7T 

3-V1+2tt 




2 <p 


> A 2 


= A 2 




lim^oo A 1/p 


1 


2 


1 + 5V2 



Table 3. The smallest eigenvalues A^ and A<s 2 under symmetry assumptions for 
the square, (a) Numerical comparison with A 2 . (b) Asymptotic behavior: the 
first row shows values of hi, the second row values of Ai for ft, f^ alf , and ftij^, 
respectively. 



A A 




Figure 10. Dependence of the numerically computed eigenvalues for the square 
on p. The three cross symbols in the graph on the right mark the values of hi for 
ft, ft^ alf , and ft£ alf . 




4.3. Rectangle. Let 

ft = {Oi,x 2 ) e R 2 \xi e (0,2), x 2 e (0,1.75)} . 

This domain was discretized using 77,312 triangles. The shape of the first eigenfunction m and the 
graph of the first eigenvalue Ai(ft;p) are similar to those for the square. However, the symmetry 
properties of the second eigenfunction 112 are different: According to our numerical observations, 
for p < 3.6 the eigenfunction 112 preserves an odd symmetry about Xi = 1 and an even symmetry 
about X2 = 0.875 (which we call Si as in the case of the square). For p > 3.7 this symmetry is 
lost and 112 maintains an odd symmetry with respect to (1,0.875), the center of ft. The contour 
lines of U2 for several values of p are shown in Fig. [l2j 
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p = 3.6 



p = 3.7 



p = 3.t 



p = S.O 



Figure 12. The numerically computed second eigenfunction ix 2 for the rectangle. 




Q = (0,2) x (0,1.75) 



p 


Ai 


A 2 


A<Si 


3.5 


12.166 


57.874 


57.874 


3.6 


12.680 


63.436 


63.436 


3.7 


13.206 


69.491 


69.492 


3.8 


13.743 


76.059 


76.083 


3.9 


14.292 


83.175 


83.256 


4.0 


14.853 


90.881 


91.059 


8.0 


49.531 


2,192.9 


2,574.6 



Figure 13. Comparison of the numerically computed second eigenvalue A 2 and 
the smallest eigenvalue Xs 1 under the symmetry Si for the rectangle ft (Ai is 
shown for reference). 



(a) 



Rectangle (0,2) x (0,1.9) 



p 


Ai 


A 2 


A<Si 


2.44 


6.5926 


20.0177 


20.0177 


2.46 


6.6579 


20.4281 


20.4281 


2.48 


6.7234 


20.8451 


20.8457 


2.50 


6.7891 


21.2683 


21.2708 


2.60 


7.1205 


23.4816 


23.5124 



(b) 



Rectangle (0,2) x (0, 1.6) 



p 


Ai 


A 2 


A<Si 


5.6 


36.077 


383.4648 


383.4648 


5.8 


38.898 


453.2332 


453.2333 


6.0 


41.892 


535.2007 


535.2009 


6.2 


45.068 


631.4438 


631.4478 


6.4 


48.437 


744.1846 


744.4026 



Table 4. Comparison of the numerically computed second eigenvalue A 2 and the 
smallest eigenvalue Xs 1 under the symmetry Si for other rectangles. 



For p > 3.7 the smallest eigenvalue Xs 1 corresponding to an eigenfunction with symmetry Si 
is larger than A 2 (cf. Fig. 13). This eigenpair can be computed by CDM on Q with an additional 
boundary condition ^(l,x 2 ) = for x 2 G (0, 1.75) or as the first eigenpair on the half-rectangle 

n haU = {(x!,x 2 ) GM 2 |xi G (0,1), x 2 G (0,1.75)}. 

Our conjecture is that for a rectangle R — (0, a) x (0, b) with < b < a there exists po > 2 such 
that U2 has two nodal domains which for p < po are rectangles with sides a/2 and b. For p > po the 
nodal domains are not rectangular and u 2 has only an odd symmetry with respect to the center 
of R. According to our numerical observations, po gets larger the larger the ratio a/b: Besides Q 
we ran the computation for two other rectangles. For R = (0,2) x (0, 1.9) the loss of symmetry 
Si of U2 is observable approximately between p = 2.44 and p = 2.48 and for R = (0, 2) x (0, 1.6) 
between p = 5.6 and p = 6.0 (cf. Table [4|. As p grows and crosses po? the nodal line which is 
straight for p < po gets distorted. This distortion is faster for smaller ratios a/b and slower for 
larger ratios. 
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4.4. Triangle with height 1. Let 

fi = {(x u x 2 ) eR 2 \ Xl e (0,l),|a? 2 | < |(1 -a?i)} 

be an isosceles triangle with base 1 and height 1. It was discretized using 38,912 triangles. Figures 
14| and [T5| show the eigenfunctions u\ and u 2 for several values of p, respectively. Table [5] lists the 



corresponding values of Ai and A 2 . 

Various intermediate path points cm were used to compute u 2 . However, the function that 
CMPA converged to did not depend on this choice. The symmetry properties of the computed u 2 
depend only on the value p. For p < 2.6 it is even in x 2 , i.e., it belongs to 

(39) <Se := {u ' — >• M | u(xi,x 2 ) = u(xi, — x 2 )}. 

For p > 2.7 this symmetry is lost by u 2 as the graphs in Fig. [T5] show. 

For p = 2.6 the computation was repeated with intermediate path points eM without symmetry 
Se but CMPA always converged to the function shown in Fig. [15] which displays symmetry Se- 

For p = 2.7 a symmetric ey[ £ <$e was chosen. The graph in Fig.[l6ja) shows how the maximum 
of the Dirichlet functional I along the path evolved during this run of CMPA. The path connecting 
u\ with —u\ which gets deformed at every step of CMPA seems to stay close to some critical point 
having symmetry Se during the first 1000 steps but then it slips down to lower values of / and 
stays close to another critical point. This is the asymmetric u 2 which the algorithm eventually 
converges to. 

Even beyond p — 2.6 there exist eigenfunctions with symmetry Se- Let u 2j s E denote a sign- 
changing eigenfunction of the p-Laplace operator on Q which lies in Se and has the smallest 
eigenvalue (which we denote A 2)4 s E ). As mentioned above, for p < 2.6 we observed that u 2 = u 2 ^s E 
(up to scaling). To compute i£2,<s E for p > 2.7 consider the following eigenvalue problem: 

-A p u = X\u\ p ~ 2 u inft half , 



(40) *£=0 onTi, 



du 

u = on T2, 



where 

ft half = {{x 1 ,x 2 )en\x2>0}, 
(4i) r 1 = {(x 1 ,x 2 )e dn haU \ x 2 = o}, 

r 2 = dn haU \ri. 



p 


Ai 


A 2 


1.1 


8.0143 


12.188 


1.2 


10.208 


16.211 


1.3 


12.673 


21.009 


1.4 


15.515 


26.847 


1.5 


18.822 


33.998 


1.6 


22.683 


42.774 


1.7 


27.196 


53.546 


1.8 


32.471 


66.762 


1.9 


38.634 


82.963 


2.0 


45.831 


102.80 


2.1 


54.228 


127.06 


2.2 


64.016 


156.72 


2.3 


75.415 


192.92 


2.4 


88.681 


237.08 


2.5 


104.10 


290.87 



P 


Ai 


A 2 




2.6 


122.02 


356.35 


356.35 


2.7 


142.81 


435.98 


435.99 


2.8 


166.94 


532.61 


532.78 


2.9 


194.90 


649.76 


650.31 


3.0 


227.29 


791.69 


792.92 


3.5 


483.05 


2,093.5 


2,107.6 


4.0 


1,006.3 


5,425.7 


5,498.4 


5.0 


4,183.4 


34,911 


35,924 


6.0 


16,688 


2.1571 • 10 5 


2.2561 • 10 5 


8.0 


2.4510 • 10 5 


7.6097 • 10 G 


8.2094 • 10 d 


10.0 


3.3583 • 10 G 


2.5069 • 10* 


2.7692 • 10* 



Table 5. Eigenvalues for the triangle with height 1. 
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p= 1.1 




p = 1.4 




p = 2.0 




p = 2.5 




p = 8.0 



Figure 14. The numerically computed first eigenfunction u\ for the triangle with 
height 1. 




p = l.l 




p = 2.6 




p = 2.7 





p = 3.0 



p = 8.0 



Figure 15. The numerically computed second eigenfunction u 2 for the triangle 
with height 1. 



(a) 



i(z max ) 



436.02 



A 2 
435.98 



u 2 ,s E p = 2.7 




(b) 

irirlr 



1000 



2000 3000 

iterations 




p = 3.5 





p = 6.0 



p= 10.0 



Figure 16. Triangle with height 1: (a) Maximum value of the Dirichlet func- 
tional I along the path during the run of CMPA for p = 2.7 and eM G <Se- (b) 
The computed eigenfunction ^2,<s E for p = 3.5, 6.0, and 10.0. 



18 



Jim horAk 



A 
100 

80 

60 

40 

20 

/ii(fi) 




1.2 
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V 




Figure 17. Dependence of the numerically computed eigenvalues for the triangle 
with height 1 on p. 



Any eigenfunction solving this problem can be extended to an eigenfunction on the whole ft by 
even symmetry about X2 — 0. Since the first eigenfunction of the original problem (pQ) for the 
triangle Q belongs to <Se, its restriction to Q haU is the first eigenfunction for (|4Q|) . Hence to 
compute ^2,<s E we just need to apply CMPA to problem (|4Q|) with paths which again connect u\ 
and —u\. The modification of the finite element method to take into account the natural boundary 
condition on Ti is straightforward. The computed values of A2,<s E are listed in Table [5] Figure 
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b) shows the corresponding eigenfunction 1^2, <s E fc> r selected values of p. 

The dependence of the eigenvalues Ai, A2, and A2,<s E on p is further plotted in Fig. [17] The 
figure also shows the limits of Ai for p — > 1 and 00 and of A2 for p — > 00 which can be computed 
explicitly. As mentioned, for example, in [18], the Cheeger constant of a triangle is given by 
fti(fi) = (Per(^) + V / 4tt|^|)/(2|^|) and hence in our case h x (ft) = 1 + v / 5 + V^tt « 5.7427. Simple 
computations yield Ai(fi) = 1 + >/5 « 3.2361 and A 2 (ft) = 1 + 9/V$ « 5.0249. 



4.5. Triangle with height 3/4. Let 

ft = {(*i,*2) G R 2 I G (0, |) , |x 2 | < § (| - a*)} 

be an isosceles triangle with base 1 and height 3/4. It was discretized using 28,672 triangles. 
Figures [18] and [19] show the eigenfunctions u\ and 112 for several values of p, respectively. Table [6] 
lists the corresponding values of Ai and A2. 

The symmetry properties of the computed 112 change again with p. For this triangle, however, 112 
gains more symmetry as p increases (unlike for the triangle with height 1 where U2 lost symmetry). 
For p < 1.6 the nodal line of U2 connects the base of the triangle with one of its other sides. For 
p > 1.7 this nodal line connects the base with the vertex above the base and 112 is odd in #2, i-e., 
it belongs to 



(42) 



Sq := {u : Q — >• R | #2) = — ^(^i, — #2)} 



p 


Ai 


A 2 


A<s 


1.1 


9.389 


14.38 


14.50 


1.2 


12.13 


19.41 


19.52 


1.3 


15.28 


25.53 


25.62 


1.4 


18.97 


33.10 


33.17 


1.5 


23.35 


42.52 


42.55 


1.6 


28.55 


54.22 


54.22 


1.7 


34.73 


68.76 


68.76 



p 


Ai 


A 2 


1.8 


42.07 


86.84 


1.9 


50.78 


109.2 


2.0 


61.11 


137.1 


2.1 


73.36 


171.6 


2.2 


87.85 


214.4 


2.3 


105.0 


267.2 


2.4 


125.2 


332.5 



p 


Ai 


A 2 


2.5 


149.1 


413.0 


3.0 


350.0 


1,196 


4.0 


1,789 


9,351 


5.0 


8,591 


6.871 • 10 4 


6.0 


3.958 • 10 4 


4.849 • 10 b 


8.0 


7.752 • 10 5 


2.232 • 10 Y 


10.0 


1.416 • 10 Y 


9.602 • 10 8 



Table 6. 



Eigenvalues for the triangle with height 3/4. 
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Figure 19. The numerically computed second eigenfunction 112 for the triangle 
with height 3/4. 




Figure 20. Higher eigenfunctions for the triangle with height 3/4 for p = 1.5: 
and ^(2,5 E ) computed as constrained local mountain pass points by CMPA with 
no a priori assumptions on symmetry, us computed by CDM enforcing symmetry 
So- 



20 



Jim horAk 



A 

120 
100 
80 
60 
40 
20 




1 




1.2 



1.4 



1.6 



1.8 



/ii(Q h 



hi(ft) 



V 




Figure 21. Dependence of the numerically computed eigenvalues for the triangle 
with height 3/4 on p. 



as can be seen in Fig. [19] It is because of lack of resolution of the numerical method close to the 
vertex (where U2 is flat) that the zero contour line in the figure for p = 1.7 and p = 8.0 does not 
exactly reach the vertex. 

For p < 1.6 there also exist eigenfunctions with symmetry So. Let Xs Q denote the smallest 
eigenvalue with an eigenfunction belonging to So (denoted us Q ). Using the notation defined in 
(jHJ) this eigenvalue can be computed using CDM on Q with an additional boundary condition 
u — on Ti or as the first eigenvalue on the half-domain Q half . The computed values of \$ are 
also listed in Table [6j For p > 1.7 the values of A2 and A<s coincide. For p = 1.6 they differ in 
the sixth digit. 

As in the previous computations, various choices of the intermediate path point eM were used 
to compute U2 on Q with no a priori assumptions on symmetry. In some cases CMPA converged 
to different functions depending on this choice (different local mountain passes). For example, for 
p = 1.5 two eigenfunctions were found: one with a nodal line connecting the base of the triangle 
with one of its sides, and another one with a nodal line connecting the two sides and having an 
even symmetry in X2- Both eigenfunctions are (numerically) local mountain pass points of / with 
respect to the constraint S. The first one is called 112 since it has the smallest eigenvalue, the 
second one is called ^(2,s E ) because of its symmetry (it could also be understood as a solution of 
(|4Q|) formulated in a similar way for the triangle with height 3/4). Both eigenfunctions are shown 



in Fig. 20 together with us Q for comparison. 



The dependence of the eigenvalues Ai, A2, and \s on p is further plotted in Fig. 21 The 



following limits of Ai and \s Q as p — > 1 and those of Ai and A2 as p — > 00 are also marked in the 
figure and have these respective values: 

fti(fi) = |(2 + v / 13 + >/6tt) ~ 6.631, Ai(fi) = §(2 + >/l3) ^ 3.737, 



ht (ft half ) = I (5 + \/T3 + 2 >/37r) ~ 9.830, A 2 (ft) = § (5 + >/l3) « 5.737. 

4.6. Equilateral triangle. For isosceles triangles close but not equal to an equilateral triangle 
a similar observation has been maded as for rectangles close but not equal to the square: the 
symmetry properties of the second eigenfunction 112 change at a certain value p ^ 2. According 
to the following computations, for an equilateral triangle this change occurs at p = 2 (as it does 
for the square). 



Let 



{(x 1 ,x 2 ) e M 2 1 x 1 e (0, ^) , \x 2 \ < 7s - ^1) } 



be an equilateral triangle with side 1. It was discretized using 32,256 triangles. With the notation 
introduced in we can define \2,s E an d ^2,<s E as in Sec. |4.4|and Xs Q and us Q as in Sec. 



Our numerical observations are summarized in Table [71 and Fig. 22 For p < 2 the second 



4.5 



eigenfunction U2 is even in X2 while for p > 2 it is odd (up to a rotation of the triangle by ±27r/3). 
We note that the values A2 listed in the table were computed with no a priori assumptions on 
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p 


Ai 


A2(= A 2 , 4 s E ) 


A<s 


1.1 


8.653 


13.37 


13.61 


1.9 


44.07 


98.20 


98.40 


2.0 


52.64 


122.8 


122.8 



P 


Ai 


A 2 (= A^) 


A2,5 E 


2.0 


52.64 


122.8 


122.8 


2.1 


62.71 


152.9 


153.2 


8.0 


4.240 • 10 5 


1.483 • 10 Y 


1.668 • 10 Y 



Table 7. Eigenvalues for the equilateral triangle with side 1. 




Figure 22. The numerically computed eigenfunctions u^, ^2,<s E anc ^ u s f° r the 
equilateral triangle. 



the symmetry of u and then compared to the computed values A2,<s E and A<s . The corresponding 
eigenfunctions 112 are in the bottom row of the figure. For p = 2 the eigenspace belonging to 
A2 is two-dimensional. The member of this eigenspace 112 to which CMPA converges depends 
on the initial path, i.e., on the choice of the intermediate point eM- The figure shows one such 
member. As already mentioned in Sec. |4.5[ it is an artifact of the numerical method that for U2 
and p = 2.1, 8.0 the zero contour line does not exactly reach the vertex, where 112 is rather flat. 



5. Remarks on the numerics 

5.1. Dependence on the mesh parameter h. Let T h denote th e set of all the triangles of a 
triangulation of Q h . The mesh parameter h was introduced in Sec. 



3.1 



as the (smallest) upper 

bound on the diameter of the circumscribed circle for triangles of T h '. In this section the depen- 
dence of the computed values of Ai and A2 on h is investigated. The investigation is conducted 
for one particular domain Q- 



-the rectangular domain used for computations in Sec. [O 

{Oi,x 2 ) e M 2 \x ± e (0,2), £ 2 e (0,1.75)} . 



Four discretizations of this domain are used. Table |8] lists details about these discretizations 
ordered by the number of triangles. Essentially, a finer mesh was obtained from a courser one by 
placing a new vertex in the middle of each triangle side of the old triangulation, in effect dividing 
each triangle in four. 

Table [9] shows the values of Ai and A2 computed for the four triangulations characterized by h 
and for selected values of p. Figure [23] gives perhaps a more telling picture: for each p it displays 
relative differences (X(T h ) — A(T' 011 ))/A(T' 011 ), where X(T h ) denotes the eigenvalue computed 
on the triangulation T h . The finest triangulation 7~' 011 is used as a reference. We observe that 
the largest differences occur for large p (here p = 8.0) and smallest differences for p close to 2. 
The differences are about twice as large for A2 computed by CMPA compared to Ai computed by 
CDM. 
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h 


number of triangles 


course 


0.079 


4,832 




0.044 


19,328 




0.022 


77,312 


fine 


0.011 


309,248 



Table 8. Triangulations used to discretize the rectangular domain Q. 



% 
1 

0.75 
0.5 
0.25 



relative difference for Ai 




% 
2 

1.5 



0.5 



relative difference for A2 

-p=l.l 
-p=1.6 

p = 4.0 / 
-p = 8.0 / 




0.011 0.022 



0.044 



0.079 





0.011 0.022 



0.044 



0.079 



Figure 23. Relative difference 



\(T h )-\(T 011 ) 



A(r .oii) 100% for Ai (left) and A 2 (right) 

computed on a triangulation T h with respect to the finest triangulation T' 011 . 
For Ai the relative differences for p = 1.1 and p = 4.0 almost coincide and cannot 
be distinguished in the graph. 





p = 1.1 


h 


Ai 


A 2 


0.079 


2.5586 


3.9507 


0.044 


2.5544 


3.9376 


0.022 


2.5533 


3.9342 


0.011 


2.5531 


3.9334 



p = 1.6 


Ai 


A 2 


4.2965 


8.2462 


4.2945 


8.2381 


4.2940 


8.2361 


4.2939 


8.2356 



p = 4.0 


Ai 


A 2 


14.884 


91.306 


14.859 


90.962 


14.853 


90.881 


14.851 


90.861 



p = S.O 


Ai 


A 2 


50.005 


2,234.8 


49.625 


2,202.7 


49.533 


2,193.6 


49.510 


2,191.0 



Table 9. Values of Ai and A2 computed for four different triangulations of the 
rectangular domain Q and p = 1.1, 1.6,4.0,8.0. 



5.2. The Augmented Lagrangian Method. As described in Sec. 3.3 this method is used to 
solve ([28]) iteratively for a given right-hand side. The Augmented Lagrangian C r defined in ([30]) 
depends on a parameter r > 0. As observed by the authors of [13J the algorithm is not very 
sensitive to the choice of r but the analysis of the influence of r on the behavior of the algorithm 
is complicated. 

The choice of r has an influence on the speed of convergence of the algorithm and at the same 
time on how precise the found numerical solutions can be. In general, for larger r the algorithm 
seems to converge faster but it is able to find only less precise approximations of the solution. For 
our computations we tried various values of r first and then chose the one which seemed to give a 
reasonable speed of convergence together with acceptable residual. This value depended strongly 
on p and also on the particular domain Q. Table [To] shows the dependence of r and of the number 
of iterations that the algorithm needed on p. For each domain one value of r was chosen from the 
given range. Similarly, the number of iterations lay in the given range. We can observe that for 
a small p a large r was needed, for a large p a smaller r. For p close to 2 we could choose r « 1. 
The number of iterations needed turned larger for p farther from 2. 
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p 


range of r 


# of iterations 


1.1 


10 4 - 10 7 


700 - 2,000 


1.2 


500 - 2, 500 


500 - 1,000 


1.8 


1 - 1.5 


80-90 


3.0 


0.3 - 0.4 


200 - 300 


10.0 


0.03 - 0.1 


1,200 - 3,000 



Table 10. The dependence of the approximate values of r and numbers of iter- 
ations in the Augmented Lagrangian Method on p. 



We note that for values of p smaller than 1.1 and larger than 10 (the particular value also de- 
pended on the domain) we were not able to find r giving satisfactory results for our implementation 
of the Augmented Lagrangian Method in conjunction with CMPA. 

5.3. CDM and CMPA. In both the Constrained Descent Method and the Constrained Mountain 
Pass Algorithm the measure of convergence is \\w u ||, the Wo' p (tt) -norm of the descent direction 
evaluated at the approximation u of the eigenfunction which is being computed. The smallest 
achieved value depended on the algorithm, on p, and in case of CMPA also on the fact whether 
there lies another critical point not far from u. For CDM the order of this value was between 10 -5 
and 10 -8 , for CMPA between 10 -3 and 10 _T . The number of iteration of CDM was approximately 
between 10 and 30. The number of iterations of CMPA varied, it depended on the shape of the 
initial path and on p, and was anywhere between 100 and 3,000. 

6. Conclusion 

In this work a concrete application of the variational numerical methods of [14] in a Banach 
space was presented. In particular, one possible choice of the descent direction required by these 
methods was proposed, implemented and tried in computations in the setting of the Sobolev space 
Wq ,p (CI). The computations yielded approximations of the smallest two Dirichlet eigenvalues and 
the corresponding eigenfunctions of the p-Laplace operator on several planar domains for p ranging 
from 1.1 to 10. This relatively large range made it possible to study the change of symmetry of 
the second eigenfunction with varying p on different domains which was first observed in (24j for 
the square and p not far from 2. The computed eigenvalues seem to agree with the asymptotic 
behavior known from theory for p — >• 1. Our range of p seems to be too small, however, in order 
to clearly observe the asymptotic behavior of the eigenpairs as p — » oo. 

Numerical experiments were conducted for the following domains: the disk, rectangles, and 
isosceles triangles. We summarize the main observations about the symmetry of the second eigen- 
function U2- For the disk it was observed that has a straight nodal line dividing the disk into 
halves for the whole range of p. 

For rectangles which are not a square and for small p the second eigenfunction is odd about its 
nodal line which is straight and connects the midpoints of the longer sides. After p crosses some 
value po > 2 there are two second eigenfunctions which are mirror images of each other. Their 
nodal line is not straight but still connects the two longer sides. 

For the square and p / 2 there are two second eigenfunctions which are images of each other 
under rotations by tt/2 about the center of the square. For p < 2 their nodal line is straight and 
connects the midpoints of the opposite sides. For p > 2 the nodal line is a diagonal of the square. 
For p = 2 there are two linearly independent second eigenfunctions. 

The symmetry observations for triangles are based on the family of isosceles triangles with 
vertices (0,-1/2), (0, 1/2), and (£, 0) with base 1 and height £ > which are symmetric about the 
xi-axis. For those shorter than the equilateral triangle and for small p there are two asymmetric 
second eigenfunctions (up to scaling) which are symmetry images of each other. Their nodal line 
connects the base with one side of the triangle. After p crosses some value po < 2 there is only 
one eigenfunction U2. It is odd about its nodal line which is straight and connects the middle of 
the base with the opposite vertex (symmetry So). 



24 



Jim horAk 



For triangles longer than the equilateral triangle and for small p there is one eigenfunction 
U2, it is even about the xi-axis (symmetry <Se) and its nodal line connects the two sides of the 
triangle. After p crosses some value po > 2 there are two asymmetric second eigenfunctions which 
are symmetry images of each other. Their nodal line still connects the two sides of the triangle. 

For the equilateral triangle and p ^ 2 there are three second eigenfunctions which are images of 
each other under rotations of the triangle about its midpoint by ±27r/3. For p < 2 their nodal line 
connects two sides of the triangle and they have even symmetry about the height coming from the 
third side. For p > 2 the nodal line of the second eigenfunctions follows a height of the triangle 
and the eigenfunctions have odd symmetry about this height. For p = 2 there are two linearly 
independent second eigenfunctions. 

Figure [20] indicates that our numerical methods could be used for finding some higher eigenfunc- 
tions and perhaps for a continuation in £ to observe the connection between these eigenfunctions 
and the second eigenfunctions for the equilateral triangle. This lies however beyond the scope of 
this paper. 

Appendix A. 

Here we give a proof of some claims used in Sections |3.2| and |3.4| A subindex notation will be 
used for general sequences and does not refer to the enumeration of eigenfunctions and eigenvalues 
in this section. 

Lemma A.l. Let (B, ||-||) be a reflexive Banach space with a strictly convex norm, I, Jg C 1 (5,]R) 
be two continuously Frechet differentiable junctionals, and u be a point in B with J(u) = 1 which 
is not a critical point of I with respect to S := {v G B | J(y) = 1}. Then the problem 

minimize L(w) := (I'(u),w) subject to w G C := {v G B | (J'(u),v) = and \\v\\ = 1} 

has a unique solution. 

Proof. First, we show that L has a negative infimum on C: L is bounded below on C by — \\I f (u)\\*. 
It attains negative values on C if there exists w G C such that L(w) ^ 0. But if L = on C, then 
we would have 

(J'(u),w)=0 => (l'(u),w) = WweB 

which would imply existence of a G R such that I'(u) — aJ'(u) = 0. This is not possible since u 
is not a critical point of I with respect to S. 

Let {w n } C C be a minimizing sequence of L, i.e., 

L(w n ) — » inf L G (—oo, 0) as n —> oo. 

Since this sequence is bounded by 1, the reflexivity of B implies existence of a subsequence (still 
denoted {w n }) which converges weakly to some w G B such that \\w\\ < 1. Since I'(u) and J'(u) 
are continuous linear functionals, we obtain 

L(w n ) —> L(w) and (J f (u),w n ) — ^ (J f (u),w) as n — » oo. 

This means that L(w) = inf^ L and (J'(u), w) = 0. 

To prove that w is a minimizer it remains to show that \\w\\ = 1. If ||u;|| < 1, then w := 
belongs to C and 

L ^ = w < L{w) 

because L(w) < 0. But this is a contradiction with the minimality of L(w). 

To show uniqueness let w\ and W2 be both minimizers. For w := \{w\ + W2) /\\\( w i + ^2)|| 
we obtain 

_ _ mine L 
w G C and minL < L(w) = -n-r- -n-. 

c ~ |||K+^2)|| 
Since the minimum is negative, this and the triangle inequality imply 

i< Uw! + \w 2 \\ < IK II + = 1. 
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Hence equality holds in the above inequalities and the strict convexity of the norm implies w\ = 
w 2 . □ 

Lemma A. 2. Let I and J be defined by (0) and S by Further, let u G S and 

(43) w u := -u + r { ,\ y- v u , where v u := (-A p ) _1 (\u\ p ~ 2 u) . 

J n \u\P~ z uv u ax v y 

TTien (I f (u), w u ) < 0. Equality holds if and only if u is a critical point of I with respect to S which 
is the case if and only if w u = 0. 

The proof of this lemma is based on the following inequality which is a direct consequence of 
the Cauchy-Schwarz and Holder inequalities. Its proof is therefore omitted. 

Auxiliary Lemma. Let f,g G Wq ,p (Q), f ^ 0. Then 

iv/l^v/v^x^n/ir 1 ^!!. 



/ 



in 

Equality holds if and only if there exists v > such that vf = g. 
Proof of Lemma \A.S\ We observe that 

(44) / \u\P- 2 uv u dx= [ (-A p v u )v u dx=\\v u \\P. 

By the definition of w u , (|44|) and the auxiliary lemma we obtain 

{I\u),w u ) = -\\u\\v + -±- [ \Vu\v~ 2 VuVv u dx 
fA*\ \\ v u\\ p Jn 

Using u G 5, testing the equation —A p v u = \u\ p ~ 2 u by u, and applying the auxiliary lemma yields 
(46) 1= / \u\"dx= [ \Vv u \ p - 2 Vv u Vudx ( *< Wvuf^Wul 



By combining ([45]) and ([46]) we conclude that {P{u), w u ) < 0. Equality holds if and only if equality 
holds in (*) and (**). According to the auxiliary lemma this is the case if and only if vu = v u for 
some v > 0. Finally, we argue that the following are equivalent: 

(a) vu = v u for some v > 0, 

(b) u is a critical point of / with respect to S, 

(c) w u 0. 

Statement (a) is equivalent to v p ~ 1 (—A p u) = \u\ p ~ 2 u and hence to (b). If (a) holds, then 
J n \u\ p ~ 2 uv u dx = v because u G S. Hence w u = — u + ^v u = and (c) holds, too. It is 
obvious that (c) implies (a). □ 



Before stating the next proposition we recall some known results (let q G (0, oo), - + ^ = 1): 



(i) The p-Laplace operator — A p : W ' p (f}) —> W~ 1,q (Q) is uniformly continuous on bounded 
sets. 

(ii) The mapping u \-> \u\ p ~ 2 u : Wq' p (Q) —> W~ 1,q (fl) is compact and uniformly continuous 
on bounded sets. 

(hi) The inverse p-Laplace operator (— A p ) _1 : W~ 1,q (Q) — >• Wq' p (Q) is uniformly continuous 
on bounded sets. 

Both claims (i) and (ii) follow from standard inequalities found, e.g., in [13, Lemmas 5.3 and 5.4]. 
The compactness in (ii) follows from the compact embedding of W ,p (tt) in L P (Q). Claim (hi) 
follows from standard inequalities found, i.e., in p~3| Propositions 5.1 and 5.2]. 
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Proposition A. 3. Let I and J be defined by (0) and S by Q). The initial value problem 

(47) j t u{t) = w u(t) , u(0) = e eS 

with w u defined in has a unique solution u(t) G S defined for t G (0, oo). There exists a 
critical point u G S of I with respect to S and a sequence {t n }^ =1 such that lim^oo t n = oo and 
lim n ^oo u(t n ) = u in Wq P (Q). 

Proof. The proof of existence of a solution and its uniqueness follows the same lines as the proof 
of Lemma 5 in |14| . Hence we focus on establishing the existence of the sequence {t n }. 

Since < I(u(T)) = I(eo) + Jq (I' \u(t)) , dt for T > and the integrand is non- 

positive, we obtain J °° \ (l'(u(i)), w u ^)\ dt < I{e$). Hence there exists a sequence {t n }^ =1 with 
lim^oo t n = oo such that for u n := u(t n ) and w n := it holds: 

(48) (7'(w n ), w n ) -» for n -» oo. 
We recall that by (g3j) and (glj) we have 

(49) w n = -u n + 1 |i p ^n 7 where v n := (-A p ) _1 (|^ n | p_2 ^n) • 

1 1 v n 1 1 



IP-2-i 



We observe that {ix n } is a bounded sequence, hence it converges weakly to some u G W ' P (Q) 
along a subsequence which we again denote {u n }. From the compactness of the map u i->> \u 
and the continuity of the inverse p-Laplacian it follows that 

(50) v n v := (-A p ) _1 (|^| p_2 ^) strongly in W^ p (Vt) as n ^ oo. 
Equation ([49]) and the fact that w n G T Un S imply 

(51) {-A p (w n + u n ),w n ) = - — 1 I \u n \ p ~ 2 u n w n dx = 0. 
Combining (|48j) and fl5T) yields 

(52) / (|V(w n + u n )\ p - 2 V(w n + u n ) - \Vu n \ p - 2 Vu n ) Vw n dx -> for ?w oo. 
On the other hand standard estimates fT3] Propositions 5.1 and 5.2] state that 

/ (\V(w n + u n )\ p - 2 V(w n + ii n )-|V^ n | p - 2 V^ n ) V^ n 

(53) ><^ ■ ""I" 112 ,, forl<p<2, 



(54) >^\\w n \\ p for2<p, 



1 

2p- 

where S > is a constant which does not depend on w n and u n . These inequalities and ([52]) imply 

(55) w n —> strongly in Wq' p (£2) as n — ^ oo. 
This and ([49]) in turn imply that {^ n } converges strongly to u and that 

(56) U =p^(-A p )- 1 (|«r 2 W ), 

which means that u is a critical point of / with respect to S. □ 

Remark. To better understand the implications of the choice of the descent direction we remark 
how the proof of the proposition would change if we used the steepest descent direction instead of 
the descent direction given by ([43]) . Up to normalization the steepest descent direction w defined 
by ([23]) can be written as the solution of 

—ApW = A p u + a\u\ p ~ 2 u 
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for a suitable a. Testing this equation by w and using w G T U S yields = — ^ {I'(u),w). 

Hence equation (|48|) would directly imply w n —> 0. We can write 

<- H^nF" 1 = || - A p Wn\\* = \\-\u n - a n \u n \ P ~ 2 U n \\^ = ^ \\l'(u n ) ~ OL n J' (u n )\\^ , 

where || • ||* denotes the norm in the dual space W~ lyq (Q). If we define H/'IsJI := inf^eR H^'C^) — 
aJ'(u)\\* as in [6], |14| . then we would obtain ||i" / |5 Un || — >• 0. The Palais-Smale condition under 
constraints which was formulated in [5] states in its simplified form that if {I'(u n )} is bounded 
and H/'Is^JI —> then {u n } possesses a convergent subsequence. In |9J it was shown that this 
condition holds in our setting. Hence the choice of the steepest descent direction would yield a 
more "classical" proof of the proposition. 
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